par(mfrow=c(2, 2))
data <- read.csv("gdp.csv")
ts_x <- ts(data$gdp)

# 时序图
plot(ts_x, type="o", pch=5, col='#39CBB4')

# 一阶差分
dif_x <- diff(ts_x)
plot(dif_x, type="o", pch=5, col='#39CBB4')

# ADF检验
# install.packages("aTSA")
library(aTSA)
adf.test(dif_x, nlag=3)

# 白噪声，纯随机检验
for( k in 1:3) print(Box.test(dif_x, lag=6*k, type="Ljung-Box"))

# 自相关图
acf(dif_x, lag.max=30)
pacf(dif_x, lag.max=30)

# 拟合 ARMA(5, 0) 模型
fit <- arima(ts_x, order=c(5, 1, 0), method="ML")


# 模型显著性检验
ts.diag(fit)

# 参数显著性检验
# 粗略的
fit
# 精确的 构造t统计量，求P值
t = abs(fit$coef)/sqrt(diag(fit$var.coef))
pt(t, length(dif_x)-length(fit$coef), lower.tail = F)

# install.packages("forecast")
library(forecast)
# 预测未来5期
forecast(fit, h = 5)


